package server.picfet;

import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import java.util.Random;

import javax.imageio.ImageIO;

import com.googlecode.javacv.cpp.opencv_core.CvMat;
import com.googlecode.javacv.cpp.opencv_core.IplImage;

public class CVisFeaExtract {
	private IplImage m_pImage;
	private int m_width;
	private int m_height;

	public double[] m_correlogram = new double[166];// 存放hsv空间中166维，距离半径为8的颜色相关图
	public double[] m_moment1L = new double[25];// 存放25个局部区域的LAB空间的L通道的一阶矩
	public double[] m_moment1A = new double[25];// 存放25个局部区域的LAB空间的A通道的一阶矩
	public double[] m_moment1B = new double[25];// 存放25个局部区域的LAB空间的B通道的一阶矩
	public double[] m_moment2L = new double[25];// 存放25个局部区域的LAB空间的L通道的二阶矩
	public double[] m_moment2A = new double[25];// 存放25个局部区域的LAB空间的A通道的二阶矩
	public double[] m_moment2B = new double[25];// 存放25个局部区域的LAB空间的B通道的二阶矩
	public double[] m_moment3L = new double[25];// 存放25个局部区域的LAB空间的L通道的三阶矩
	public double[] m_moment3A = new double[25];// 存放25个局部区域的LAB空间的A通道的三阶矩
	public double[] m_moment3B = new double[25];// 存放25个局部区域的LAB空间的B通道的三阶矩
	public int[] m_hist = new int[166];// 存放HSV空间的166维的颜色直方图
	public int[][][] m_EdgeHistogram = new int[5][8][8];// 存放5个区域的量化成8×8的边缘直方图
	public double[] m_Energy = new double[24];// 存放24个方向的共生矩阵上提取的energy值
	public double[] m_Entropy = new double[24];// 存放24个方向的共生矩阵上提取的entropy值
	public double[] m_Contrast = new double[24];// 存放24个方向的共生矩阵上提取的contrast值
	public double[] m_Correlation = new double[24];// 存放24个方向的共生矩阵上提取的correlation值
	public double[] m_waveletfeature = new double[108];// 存放3×3分区进行harr小波变换得到的108个纹理特征

	public CVisFeaExtract() {
		m_pImage = null;
		m_width = 0;
		m_height = 0;
	}

	public CVisFeaExtract(IplImage pImage) {
		if (pImage != null) {
			m_pImage = pImage;
			m_width = pImage.width();
			m_height = pImage.height();
		} else {
			m_pImage = null;
			m_width = 0;
			m_height = 0;
		}
	}

	// 加载一张图片
	public boolean LoadImage(String fileName) throws IOException {
		File file = new File(fileName);
		if (!file.isFile()) {
			return false;
		}
		BufferedImage img = ImageIO.read(file);
		IplImage origImg = IplImage.createFrom(img);
		if (origImg == null)
			return false;
		this.m_pImage = origImg;
		this.m_width = origImg.width();
		this.m_height = origImg.height();
		return true;
	}

	public void ReleaseImage() {
		if (this.m_pImage != null) {
			m_pImage.release();
			this.m_width = 0;
			this.m_height = 0;
		}
	}

	// 返回图片的类型
	public int GetTypeFromFileName(String fileName) {
		int type = 0;

		return type;
	}

	
	public boolean ExtractColorCorrelogram() {
		CvMat mat = m_pImage.asCvMat();
		int colorarray[][] = new int[this.m_width][this.m_height];
		int color;
		int hist[] = new int[166];// 临时存放颜色直方图
		double correlogram[][] = new double[166][8];

		// 对存放颜色相关图的数组correlogram等进行初始化(该相关图是颜色自相关性，d＝8）
		for (int i = 0; i < m_width; i++)
			for (int j = 0; j < m_height; j++) {
				colorarray[i][j] = 0;
			}

		for (int i = 0; i < 166; i++)
			for (int j = 0; j < 8; j++)
				correlogram[i][j] = 0.0;

		for (int i = 0; i < 166; i++)
			m_correlogram[i] = 0.0;

		for (int i = 0; i < 166; i++)
			hist[i] = 0;

		for (int i = 0; i < m_width; i++)
			for (int j = 0; j < m_height; j++) {
				// 取一个象素的RGB值

				int r = (int) mat.get(j, i, 0);
				int g = (int) mat.get(j, i, 1);
				int b = (int) mat.get(j, i, 2);

				// 将RGB转化为166维的HSV
				color = RGB2HSV166(r, g, b);

				// 最后输出166维的直方图
				hist[color]++;

				colorarray[i][j] = color;
			}

		for (int k = 1; k <= 8; k++) {
			for (int i = 0; i < m_width; i++)
				for (int j = 0; j < m_height; j++) {
					color = colorarray[i][j];
					correlogram[color][k - 1] += func(i - k, j - k + 1,
							2 * k - 2, 0, color, colorarray)
							+ func(i - k, j - k, 2 * k, 1, color, colorarray)
							+ func(i - k, j + k, 2 * k, 1, color, colorarray)
							+ func(i + k, j - k + 1, 2 * k - 2, 0, color,
									colorarray);
				}
		}

		for (int i = 0; i < 166; i++) {
			if (hist[i] <= 0)
				hist[i] = 1;
			for (int k = 1; k <= 8; k++) {
				correlogram[i][k - 1] = correlogram[i][k - 1]
						/ (double) (hist[i] * 8 * k);
				m_correlogram[i] += correlogram[i][k - 1];
			}
		}
		return false;
	}

	private int RGB2HSV166(double red, double green, double blue) {
		// RGB转化为HSV

		double h, s, v;
		int hlevel, slevel = 0, vlevel = 0;
		double max, min;
		int result;
		double r, g, b;
		r = red / 255.0;
		g = green / 255.0;
		b = blue / 255.0;
		max = max(max(r, g), max(g, b));
		min = min(min(r, g), min(g, b));
		v = max;

		if (max == 0) {
			s = 0;
			h = 0;
		} else {
			s = (max - min) / max;
			if (max != min) {
				if (r == max) {
					h = (g - b) / (max - min);
				} else {
					if (g == max) {
						h = 2.0 + (b - r) / (max - min);
					} else {
						h = 4.0 + (r - g) / (max - min);
					}
				}
				h = h * 60.0;
				if (h < 0)
					h = h + 360.0;
			} else
				h = 0;

		}

		// 将HSV空间量化为166为颜色向量（18个h×3个s×3个v＋4个gray）
		int gray = 0;
		int f;

		// 4个灰度极
		if (s < 0.1 && v > 0.1 && v <= 0.4)
			gray = 1;
		else if (s < 0.1 && v > 0.4 && v <= 0.7)
			gray = 2;
		else if (s < 0.1 && v > 0.7 && v <= 1.0)
			gray = 3;

		hlevel = (int) h / 20;

		if ((s >= 0.1) && (s <= 0.4))
			slevel = 0;
		else if ((s > 0.4) && (s <= 0.7))
			slevel = 1;
		else if ((s > 0.7) && (s <= 1))
			slevel = 2;

		if ((v >= 0.1) && (v <= 0.4))
			vlevel = 0;
		else if ((v > 0.4) && (v <= 0.7))
			vlevel = 1;
		else if ((v > 0.7) && (v <= 1))
			vlevel = 2;

		// 输出量化成166维的颜色
		if (v <= 0.1)
			result = 0;
		else if (s < 0.1 && v > 0.1)
			result = gray;
		else {
			f = 3 * 3 * hlevel + 3 * slevel + vlevel + 4;// f的值为4－165
			result = f;
		}
		return result;
	}

	// 一个递归函数，统计与点（x，y）距离为k的象素个数，flag＝1为水平方向统计；flag＝0为垂直方向统计
	int func(int x, int y, int k, int flag, int color, int[][] colorarray) {
		// 判断是否越界
		if (x >= m_width || y >= m_height || x < 0 || y < 0)
			return 0;

		if (k == 0) {
			if (colorarray[x][y] == color)
				return 1;
			else
				return 0;
		} else {
			if (flag == 1)
				return func(x, y, k - 1, flag, color, colorarray)
						+ func(x + k, y, 0, flag, color, colorarray);
			else if (flag == 0)
				return func(x, y, k - 1, flag, color, colorarray)
						+ func(x, y + k, 0, flag, color, colorarray);
		}
		return 0;
	}

	public boolean ExtractColorMoments() {
		CvMat mat = m_pImage.asCvMat();
		double colorarray[][][] = new double[m_width][m_height][3];
		double Lvalue;
		double Avalue;
		double Bvalue;

		// 把每个象素的LAB值存放在colorarray数组里
		for (int i = 0; i < m_width; i++)
			for (int j = 0; j < m_height; j++) {
				// 取一个象素的RGB值
				int r = (int) mat.get(j, i, 0);
				int g = (int) mat.get(j, i, 1);
				int b = (int) mat.get(j, i, 2);

				colorarray[i][j][0] = RGB2LAB(r, g, b, 1);
				colorarray[i][j][1] = RGB2LAB(r, g, b, 2);
				colorarray[i][j][2] = RGB2LAB(r, g, b, 3);
			}

		// 分块计算一阶矩
		for (int m = 0; m < 5; m++)
			for (int n = 0; n < 5; n++) {
				Lvalue = 0.0;
				Avalue = 0.0;
				Bvalue = 0.0;
				int N = 0;
				for (int i = m * (m_height / 5); i < (m + 1) * (m_height / 5)
						&& i < m_height; i++)
					for (int j = n * (m_width / 5); j < (n + 1) * (m_width / 5)
							&& j < m_width; j++) {
						Lvalue += colorarray[j][i][0];
						Avalue += colorarray[j][i][1];
						Bvalue += colorarray[j][i][2];
						N++;
					}

				m_moment1L[m * 5 + n] = Lvalue / N;
				m_moment1A[m * 5 + n] = Avalue / N;
				m_moment1B[m * 5 + n] = Bvalue / N;
			}

		// 分块计算二阶矩
		for (int m = 0; m < 5; m++)
			for (int n = 0; n < 5; n++) {
				Lvalue = 0.0;
				Avalue = 0.0;
				Bvalue = 0.0;
				int N = 0;
				for (int i = m * (m_height / 5); i < (m + 1) * (m_height / 5)
						&& i < m_height; i++)
					for (int j = n * (m_width / 5); j < (n + 1) * (m_width / 5)
							&& j < m_width; j++) {
						Lvalue += Math.pow((colorarray[j][i][0] - m_moment1L[m
								* 5 + n]), 2.0);
						Avalue += Math.pow((colorarray[j][i][1] - m_moment1A[m
								* 5 + n]), 2.0);
						Bvalue += Math.pow((colorarray[j][i][2] - m_moment1B[m
								* 5 + n]), 2.0);
						N++;
					}
				m_moment2L[m * 5 + n] = Math.sqrt(Lvalue / N);
				m_moment2A[m * 5 + n] = Math.sqrt(Avalue / N);
				m_moment2B[m * 5 + n] = Math.sqrt(Bvalue / N);
			}

		// 分块计算三阶矩
		for (int m = 0; m < 5; m++)
			for (int n = 0; n < 5; n++) {
				Lvalue = 0.0;
				Avalue = 0.0;
				Bvalue = 0.0;
				int N = 0;
				for (int i = m * (m_height / 5); i < (m + 1) * (m_height / 5)
						&& i < m_height; i++)
					for (int j = n * (m_width / 5); j < (n + 1) * (m_width / 5)
							&& j < m_width; j++) {
						Lvalue += Math.pow((colorarray[j][i][0] - m_moment1L[m
								* 5 + n]), 3.0);
						Avalue += Math.pow((colorarray[j][i][1] - m_moment1A[m
								* 5 + n]), 3.0);
						Bvalue += Math.pow((colorarray[j][i][2] - m_moment1B[m
								* 5 + n]), 3.0);
						N++;
					}

				m_moment3L[m * 5 + n] = Math.pow(Math.abs(Lvalue / N),
						1.0 / 3.0);
				m_moment3A[m * 5 + n] = Math.pow(Math.abs(Avalue / N),
						1.0 / 3.0);
				m_moment3B[m * 5 + n] = Math.pow(Math.abs(Bvalue / N),
						1.0 / 3.0);
				// 这是为了解决数学函数异常问题：pow函数的第一参数为负数，第二参数为小数时输出异常
				// caojuan2006-09-12
				if (Lvalue < 0)
					m_moment3L[m * 5 + n] = -m_moment3L[m * 5 + n];
				if (Avalue < 0)
					m_moment3A[m * 5 + n] = -m_moment3A[m * 5 + n];
				if (Bvalue < 0)
					m_moment3B[m * 5 + n] = -m_moment3B[m * 5 + n];

			}
		return true;
	}

	// flag = 1 return L,flag = 2 return A,flag = 3 return B.
	private double RGB2LAB(int rn, int gn, int bn, int flag) {
		double x, y, z;
		double rm = rn / 255.0;
		double gm = gn / 255.0;
		double bm = bn / 255.0;
		double L, A, B;
		x = 0.4124 * rm + 0.3576 * gm + 0.1805 * bm;
		y = 0.2126 * rm + 0.7152 * gm + 0.0722 * bm;
		z = 0.0193 * rm + 0.1192 * gm + 0.9505 * bm;

		// 对于标准的CIE光照D65，参考白色x＝95.05;y=100;z=108.91
		x = x / 0.9505;
		y = y;
		z = z / 1.0891;
		if (y < 0.008856)
			L = 903.3 * y;
		else
			L = 116 * Math.pow(y, 1.0 / 3) - 16;
		A = 500 * (f(x) - f(y));
		B = 200 * (f(y) - f(z));

		return flag == 1 ? L : (flag == 2 ? A : (flag == 3 ? B : 0));
	}

	double f(double w) {
		double result;
		if (w > 0.008856)
			result = Math.pow(w, 1.0 / 3);
		else
			result = 7.787 * w + 0.137931;
		return result;
	}

	// *******************************************************************************************//
	// Computing ColorHistogram
	// *******************************************************************************************//
	// 该函数输入一帧图像m_image，然后输出它在HSV空间上的166维直方图hist[]
	public boolean ExtractHistogram() {

		int color;
		CvMat mat = m_pImage.asCvMat();
		// 给直方图数组hist初始化
		//!此处原来所提供的代码错误 数组越界 for (int i = 0; i <= 166; i++)
		for (int i = 0; i < 166; i++)
			m_hist[i] = 0;

		for (int i = 0; i < mat.rows(); i++)
			for (int j = 0; j < mat.cols(); j++) {
				// 取一个象素的RGB值

				int b = (int) mat.get(i, j, 0);
				int g = (int) mat.get(i, j, 1);
				int r = (int) mat.get(i, j, 2);

				color = RGB2HSV166(r, g, b);

				// 最后输出166维的直方图
				m_hist[color]++;

			}

		return true;
	}

	// *******************************************************************************************//
	// Computing EdgeHistogram
	// *******************************************************************************************//
	// 用sobel算子提取边缘直方图；分成五块提取，四个角，以及中间的一个重叠块

	public boolean ExtractEdgeHistogram() {
		int mask0[] = { -1, -2, -1, 0, 0, 0, 1, 2, 1 };// 0度
		int mask1[] = { -2, -1, 0, -1, 0, 1, 0, 1, 2 };// 45度
		int mask2[] = { -1, 0, 1, -2, 0, 2, -1, 0, 1 };// 90度
		int mask3[] = { 0, 1, 2, -1, 0, 1, -2, -1, 0 };// 135度
		int mask4[] = { 1, 2, 1, 0, 0, 0, -1, -2, -1 };// 180度
		int mask5[] = { 2, 1, 0, 1, 0, -1, 0, -1, -2 };// 225度
		int mask6[] = { 1, 0, -1, 2, 0, -2, 1, 0, -1 };// 270度
		int mask7[] = { 0, -1, -2, 1, 0, -1, 2, 1, 0 };// 315度
		double result[] = new double[8];// 存放sobel算子提取边缘的中间结果

		double gray;
		CvMat mat = m_pImage.asCvMat();
		// 初始化直方图数组
		for (int d = 0; d < 5; d++)
			for (int i = 0; i < 8; i++)
				for (int j = 0; j < 8; j++)
					m_EdgeHistogram[d][i][j] = 0;

		// 计算r0区域的八个方向的边缘

		for (int x = m_height / 4; x < 3 * m_height / 4; x++)
			for (int y = m_width / 4; y < 3 * m_width / 4; y++) {
				for (int k = 0; k < 8; k++)
					result[k] = 0;

				for (int j = -1; j <= 1; j++) {
					for (int i = -1; i <= 1; i++) {
						// gray=(double)m_pImage->GetPixelGray(y+i,x+j);
//						int r = (int) mat.get(i + y, j + x, 0);
//						int g = (int) mat.get(i + y, j + x, 1);
//						int b = (int) mat.get(i + y, j + x, 2);
						int r = (int) mat.get(j + x, i + y, 0);
						int g = (int) mat.get(j + x, i + y, 1);
						int b = (int) mat.get(j + x, i + y, 2);
						// 计算灰度值
						gray = r * 0.3 + g * 0.59 + b * 0.11;
						result[0] += mask0[(j + 1) * 3 + i + 1] * gray;
						result[1] += mask1[(j + 1) * 3 + i + 1] * gray;
						result[2] += mask2[(j + 1) * 3 + i + 1] * gray;
						result[3] += mask3[(j + 1) * 3 + i + 1] * gray;
						result[4] += mask4[(j + 1) * 3 + i + 1] * gray;
						result[5] += mask5[(j + 1) * 3 + i + 1] * gray;
						result[6] += mask6[(j + 1) * 3 + i + 1] * gray;
						result[7] += mask7[(j + 1) * 3 + i + 1] * gray;
					}
				}
				// 取值最大的方向作为该象素点的方向
				int maxresult = (int) result[0];
				int maxk = 0;
				for (int k = 1; k < 8; k++) {
					if (result[k] > maxresult)

					{
						maxresult = (int) result[k];
						maxk = k;
					}
				}
				if (maxresult < 0)
					maxresult = 0;
				maxresult = maxresult * 8 / (256 * 4);// 将最大方向上的sobel结果量化成8级
				m_EdgeHistogram[0][maxk][maxresult]++;
			}

		// 计算r1区域的八个方向的边缘

		for (int x = 1; x < m_height / 2; x++)
			for (int y = 1; y < m_width / 2; y++) {
				for (int k = 0; k < 8; k++)
					result[k] = 0;
				for (int j = -1; j <= 1; j++) {
					for (int i = -1; i <= 1; i++) {
						// gray=m_pImage->GetPixelGray(y+i,x+j);
//						int r = (int) mat.get(i + y, j + x, 0);
//						int g = (int) mat.get(i + y, j + x, 1);
//						int b = (int) mat.get(i + y, j + x, 2);
						int r = (int) mat.get(j + x, i + y, 0);
						int g = (int) mat.get(j + x, i + y, 1);
						int b = (int) mat.get(j + x, i + y, 2);
						// 计算灰度值
						gray = r * 0.3 + g * 0.59 + b * 0.11;
						result[0] += mask0[(j + 1) * 3 + i + 1] * gray;
						result[1] += mask1[(j + 1) * 3 + i + 1] * gray;
						result[2] += mask2[(j + 1) * 3 + i + 1] * gray;
						result[3] += mask3[(j + 1) * 3 + i + 1] * gray;
						result[4] += mask4[(j + 1) * 3 + i + 1] * gray;
						result[5] += mask5[(j + 1) * 3 + i + 1] * gray;
						result[6] += mask6[(j + 1) * 3 + i + 1] * gray;
						result[7] += mask7[(j + 1) * 3 + i + 1] * gray;
					}
				}
				// 取值最大的方向作为该象素点的方向
				int maxresult = (int) result[0];
				int maxk = 0;
				for (int k = 1; k < 8; k++) {
					if (result[k] > maxresult)

					{
						maxresult = (int) result[k];
						maxk = k;
					}
				}
				if (maxresult < 0)
					maxresult = 0;
				maxresult = maxresult * 8 / (256 * 4);// 将最大方向上的sobel结果量化成8级
				m_EdgeHistogram[1][maxk][maxresult]++;
			}

		// 计算r2区域的八个方向的边缘

		for (int x = 1; x < m_height / 2; x++)
			for (int y = m_width / 2; y < m_width - 1; y++) {
				for (int k = 0; k < 8; k++)
					result[k] = 0;

				for (int j = -1; j <= 1; j++) {
					for (int i = -1; i <= 1; i++) {
						// gray=m_pImage->GetPixelGray(y+i,x+j);
//						int r = (int) mat.get(i + y, j + x, 0);
//						int g = (int) mat.get(i + y, j + x, 1);
//						int b = (int) mat.get(i + y, j + x, 2);
						int r = (int) mat.get(j + x, i + y, 0);
						int g = (int) mat.get(j + x, i + y, 1);
						int b = (int) mat.get(j + x, i + y, 2);
						// 计算灰度值
						gray = r * 0.3 + g * 0.59 + b * 0.11;
						result[0] += mask0[(j + 1) * 3 + i + 1] * gray;
						result[1] += mask1[(j + 1) * 3 + i + 1] * gray;
						result[2] += mask2[(j + 1) * 3 + i + 1] * gray;
						result[3] += mask3[(j + 1) * 3 + i + 1] * gray;
						result[4] += mask4[(j + 1) * 3 + i + 1] * gray;
						result[5] += mask5[(j + 1) * 3 + i + 1] * gray;
						result[6] += mask6[(j + 1) * 3 + i + 1] * gray;
						result[7] += mask7[(j + 1) * 3 + i + 1] * gray;
					}
				}
				// 取值最大的方向作为该象素点的方向
				int maxresult = (int) result[0];
				int maxk = 0;
				for (int k = 1; k < 8; k++) {
					if (result[k] > maxresult)

					{
						maxresult = (int) result[k];
						maxk = k;
					}
				}
				if (maxresult < 0)
					maxresult = 0;
				maxresult = maxresult * 8 / (256 * 4);// 将最大方向上的sobel结果量化成8级
				m_EdgeHistogram[2][maxk][maxresult]++;

			}

		// 计算r3区域的八个方向的边缘

		for (int x = m_height / 2; x < m_height - 1; x++)
			for (int y = 1; y < m_width / 2; y++) {
				for (int k = 0; k < 8; k++)
					result[k] = 0;

				for (int j = -1; j <= 1; j++) {
					for (int i = -1; i <= 1; i++) {
						// gray=m_pImage->GetPixelGray(y+i,x+j);
						int r = (int) mat.get(i + y, j + x, 0);
						int g = (int) mat.get(i + y, j + x, 1);
						int b = (int) mat.get(i + y, j + x, 2);
						// 计算灰度值
						gray = r * 0.3 + g * 0.59 + b * 0.11;
						result[0] += mask0[(j + 1) * 3 + i + 1] * gray;
						result[1] += mask1[(j + 1) * 3 + i + 1] * gray;
						result[2] += mask2[(j + 1) * 3 + i + 1] * gray;
						result[3] += mask3[(j + 1) * 3 + i + 1] * gray;
						result[4] += mask4[(j + 1) * 3 + i + 1] * gray;
						result[5] += mask5[(j + 1) * 3 + i + 1] * gray;
						result[6] += mask6[(j + 1) * 3 + i + 1] * gray;
						result[7] += mask7[(j + 1) * 3 + i + 1] * gray;
					}
				}
				// 取值最大的方向作为该象素点的方向
				int maxresult = (int) result[0];
				int maxk = 0;
				for (int k = 1; k < 8; k++) {
					if (result[k] > maxresult)

					{
						maxresult = (int) result[k];
						maxk = k;
					}
				}
				if (maxresult < 0)
					maxresult = 0;
				maxresult = maxresult * 8 / (256 * 4);// 将最大方向上的sobel结果量化成8级
				m_EdgeHistogram[3][maxk][maxresult]++;

			}

		// 计算r4区域的八个方向的边缘

		for (int x = m_height / 2; x < m_height - 1; x++)
			for (int y = m_width / 2; y < m_width - 1; y++) {
				for (int k = 0; k < 8; k++)
					result[k] = 0;

				for (int j = -1; j <= 1; j++) {
					for (int i = -1; i <= 1; i++) {
						// gray=m_pImage->GetPixelGray(y+i,x+j);
//						int r = (int) mat.get(i + y, j + x, 0);
//						int g = (int) mat.get(i + y, j + x, 1);
//						int b = (int) mat.get(i + y, j + x, 2);
						int r = (int) mat.get(j + x, i + y, 0);
						int g = (int) mat.get(j + x, i + y, 1);
						int b = (int) mat.get(j + x, i + y, 2);
						// 计算灰度值
						gray = r * 0.3 + g * 0.59 + b * 0.11;
						result[0] += mask0[(j + 1) * 3 + i + 1] * gray;
						result[1] += mask1[(j + 1) * 3 + i + 1] * gray;
						result[2] += mask2[(j + 1) * 3 + i + 1] * gray;
						result[3] += mask3[(j + 1) * 3 + i + 1] * gray;
						result[4] += mask4[(j + 1) * 3 + i + 1] * gray;
						result[5] += mask5[(j + 1) * 3 + i + 1] * gray;
						result[6] += mask6[(j + 1) * 3 + i + 1] * gray;
						result[7] += mask7[(j + 1) * 3 + i + 1] * gray;
					}
				}
				// 取值最大的方向作为该象素点的方向
				int maxresult = (int) result[0];
				int maxk = 0;
				for (int k = 1; k < 8; k++) {
					if (result[k] > maxresult)

					{
						maxresult = (int) result[k];
						maxk = k;
					}
				}
				if (maxresult < 0)
					maxresult = 0;
				maxresult = maxresult * 8 / (256 * 4);// 将最大方向上的sobel结果量化成8级
				m_EdgeHistogram[4][maxk][maxresult]++;

			}
		return false;
	}

	// *******************************************************************************************//
	// Computing TextureCooccurence
	// *******************************************************************************************//
	// 从24个灰度级的共生矩阵中提取四个特征：entropy，energy，contrast，correlation；分别将结果存放在
	// m_Entropy[24],m_Energy[24];m_Contrast[24];m_Correlation[24]中
	public boolean ExtractTextureCooccurence() {
		double[][][] PMatrix = new double[24][256][256];
		for (int k = 0; k < 24; k++) {
			for (int i = 0; i < 256; i++) {
				for (int j = 0; j < 256; j++) {
					PMatrix[k][i][j] = 0.0;
				}
			}
		}

		ComputeMatrix(PMatrix);

		// 对共生矩阵PMatrix进行归一化

		for (int k = 0; k < 24; k++) {
			double total = 0.0;
			for (int i = 0; i < 256; i++)
				for (int j = 0; j < 256; j++)
					total += PMatrix[k][i][j];
			for (int i = 0; i < 256; i++)
				for (int j = 0; j < 256; j++) {
					if (total != 0)
						PMatrix[k][i][j] = (double) PMatrix[k][i][j]
								/ (double) total;
					else
						PMatrix[k][i][j] = 0.0;
				}
		}

		for (int i = 0; i < 24; i++) {
			m_Energy[i] = 0.0;
			m_Entropy[i] = 0.0;
			m_Contrast[i] = 0.0;
			m_Correlation[i] = 0.0;
		}

		// 计算24个共生矩阵上的四个特征
		for (int k = 0; k < 24; k++) {
			for (int i = 0; i < 256; i++)
				for (int j = 0; j < 256; j++) {
					// 能量
					m_Energy[k] += PMatrix[k][i][j] * PMatrix[k][i][j];

					// 熵
					if (PMatrix[k][i][j] > 1e-12) {
						m_Entropy[k] -= PMatrix[k][i][j]
								* Math.log((double) PMatrix[k][i][j]);
					}

					// contrast
					m_Contrast[k] += (double) (i - j) * (double) (i - j)
							* PMatrix[k][i][j];

				}

			// 计算ux
			double ux = 0.0;
			double localtotal = 0.0;
			for (int i = 0; i < 256; i++) {
				localtotal = 0.0;
				for (int j = 0; j < 256; j++) {
					localtotal += (double) PMatrix[k][i][j];
				}
				ux += (double) i * localtotal;
			}

			// 计算uy
			double uy = 0.0;
			for (int j = 0; j < 256; j++) {
				localtotal = 0.0;
				for (int i = 0; i < 256; i++) {
					localtotal += (double) PMatrix[k][i][j];
				}
				uy += (double) j * localtotal;
			}

			// 计算sigmax
			double sigmax = 0.0;
			for (int i = 0; i < 256; i++) {
				localtotal = 0.0;
				for (int j = 0; j < 256; j++) {
					localtotal += (double) PMatrix[k][i][j];
				}
				sigmax += (double) (i - ux) * (double) (i - ux) * localtotal;
			}

			// 计算sigmay
			double sigmay = 0.0;
			for (int j = 0; j < 256; j++) {
				localtotal = 0.0;
				for (int i = 0; i < 256; i++) {
					localtotal += (double) PMatrix[k][i][j];
				}
				sigmay += (double) (j - uy) * (double) (j - uy) * localtotal;
			}

			// 计算相关性
			m_Correlation[k] = 0.0;
			for (int i = 0; i < 256; i++) {
				for (int j = 0; j < 256; j++) {
					m_Correlation[k] += (double) (i - ux) * (double) (j - uy)
							* PMatrix[k][i][j];
				}
			}
			if (sigmax == 0)
				m_Correlation[k] = 0;
			else
				m_Correlation[k] /= sigmax;
			m_Correlation[k] /= sigmay;
		}

		return true;

	}

	public boolean ComputeMatrix(double[][][] PMatrix) {
		int[][] image = new int[m_width][m_height];
		CvMat mat = m_pImage.asCvMat();
		// 将彩色图像量化为256个灰度级
		for (int i = 0; i < m_width; i++)
			for (int j = 0; j < m_height; j++) {
				// image[i][j]=m_pImage->GetPixelGray(i,j);
				int r = (int) mat.get(i, j, 0);
				int g = (int) mat.get(i, j, 1);
				int b = (int) mat.get(i, j, 2);
				// 计算灰度值
				image[i][j] = (int) (r * 0.3 + g * 0.59 + b * 0.11);

			}

		// 计算90度的距离分别从1到6的6个灰度共现阵
		for (int d = 1; d <= 6; d++) {
			for (int i = 0; i < m_width; i++) {
				for (int j = 0; j < m_height - d; j++) {
					PMatrix[2 * 6 + d - 1][image[i][j]][image[i][j + d]] += 1;
					PMatrix[2 * 6 + d - 1][image[i][j + d]][image[i][j]] += 1;
				}
			}
		}

		// 计算0度的灰度共现阵
		for (int d = 1; d <= 6; d++) {
			for (int i = 0; i < m_width - d; i++) {
				for (int j = 0; j < m_height; j++) {
					PMatrix[d - 1][image[i][j]][image[i + d][j]] += 1;
					PMatrix[d - 1][image[i + d][j]][image[i][j]] += 1;
				}
			}
		}

		// 计算45度的灰度共现阵
		for (int d = 1; d <= 6; d++) {
			for (int i = 0; i < m_width - d; i++) {
				for (int j = 0; j < m_height - d; j++) {
					PMatrix[6 + d - 1][image[i][j]][image[i + d][j + d]] += 1;
					PMatrix[6 + d - 1][image[i + d][j + d]][image[i][j]] += 1;
				}
			}
		}

		// 计算135度的灰度共现阵
		for (int d = 1; d <= 6; d++) {
			for (int i = d; i < m_width; i++) {
				for (int j = 0; j < m_height - d; j++) {
					PMatrix[3 * 6 + d - 1][image[i][j]][image[i - d][j + d]] += 1;
					PMatrix[3 * 6 + d - 1][image[i - d][j + d]][image[i][j]] += 1;
				}
			}
		}

		return true;
	}

	// *******************************************************************************************//
	// Computing HaarWavelet
	// *******************************************************************************************//
	public boolean ExtractHaarWavelet() {
		short[][] spOriginData = new short[m_height][m_width];
		short[][] spTransData0 = new short[m_height][m_width];
		short[][] spTransData1 = new short[m_height][m_width];
		CvMat mat = m_pImage.asCvMat();
		// 将原始图像变为灰度图
		for (int i = 0; i < m_height; i++)
			for (int j = 0; j < m_width; j++) {
				int r = (int) mat.get(i, j, 0);
				int g = (int) mat.get(i, j, 1);
				int b = (int) mat.get(i, j, 2);
				// 计算灰度值
				spOriginData[i][j] = (short) (0.3 * r + 0.59 * g + 0.11 * b);
			}

		// 将原始图像划分为3×3的格子，分别对每个格子进行三层haar小波变换
		int intervalH = m_height / 3;
		int intervalW = m_width / 3;
		for (int i = 0; i < 3; i++)
			for (int j = 0; j < 3; j++) {
				// (hstart,wstart)为被处理区域的原点
				int Hstart = intervalH * i;
				int Wstart = intervalW * j;
				double Vr[] = new double[12];
				Haar_MultiLayers(spOriginData, spTransData0, spTransData1,
						Hstart, Wstart, intervalH, intervalH / 2, intervalW,
						intervalW / 2, 3, Vr);
				for (int k = 0; k < 12; k++)
					m_waveletfeature[(i * 3 + j) * 12 + k] = Vr[k];
			}
		return true;
	}

	/********************************************************************************
	 * 函数描述： Haar_MultiLayers完成多层图像的小波变换 * 函数参数： short
	 * **spOriginData：二维指针，指向原始的图像数据 * short
	 * **spTransData0：小波变换系数，存放一次水平变换后的小波系数 * short
	 * **spTransData1：小波变换系数，存放一次数值变换后的小波系数 * int Hstart ：图像属性参数，被处理图像的y轴起始位置 *
	 * int Wstart ：图像属性参数，被处理图像的x轴起始位置 * int nHeight ：图像属性参数，数值为原始图像的高度值 * int
	 * nHeight_H ：图像属性参数，数值为原始图像高度值的一半 * int nWidth ：图像属性参数，数值为原始图像的宽度值 * int
	 * nWidth_H ：图像属性参数，数值为原始图像宽度值的一半 * int layer ：小波变换的层数, * double Vr
	 * ：返回三层小波变换后，每个子频带的方差 *
	 ********************************************************************************/
	public boolean Haar_MultiLayers(short[][] spOriginData,
			short[][] spTransData0, short[][] spTransData1, int Hstart,
			int Wstart, int nHeight, int nHeight_H, int nWidth, int nWidth_H,
			int layer, double[] Vr) {

		int i;
		short[][] pData, pTran0, pTran1;
		// 图像的属性参数
		int iWidth, iHeight, iWidth_H, iHeight_H, iHstart, iWstart;
		// 获得数据空间的指针
		pData = spOriginData;
		pTran0 = spTransData0;
		pTran1 = spTransData1;
		// 图像属性参数赋值
		iHstart = Hstart;
		iWstart = Wstart;
		iWidth = nWidth;
		iWidth_H = nWidth_H;
		iHeight = nHeight;
		iHeight_H = nHeight_H;
		// 利用循环完成三次小波变换
		for (i = 1; i <= layer; i++) {
			Haar_Once(pData, pTran0, pTran1, iHstart, iWstart, iHeight,
					iHeight_H, iWidth, iWidth_H, i);
			// 由于小波系数可能不在0－255的范围内，因此对一次haar小波变换后的结果pTran1映射到0－255，变暂时保存在spData中
			// 供计算variance时用
			// ***********************将小波系数影射到0-255的可显示范围*********************
			// 计算得到图像小波系数的极大值与极小值
			short MaxPixVal, MinPixVal, Diff;
			MaxPixVal = pTran1[0][0];
			MinPixVal = pTran1[0][0];
			for (int y = 0; y < iHeight; y++) {
				for (int x = 0; x < iWidth; x++) {
					if (MaxPixVal < pTran1[iHstart + y][iWstart + x])
						MaxPixVal = pTran1[iHstart + y][iWstart + x];
					if (MinPixVal > pTran1[iHstart + y][iWstart + x])
						MinPixVal = pTran1[iHstart + y][iWstart + x];
				}
			}
			// 计算获得小波系数的极值差
			Diff = (short) (MaxPixVal - MinPixVal);

			for (int y = 0; y < iHeight; y++)
				for (int x = 0; x < iWidth; x++)

				{// 因为小波变换后的小波系数有可能超过255甚至更多，那么就将
					// 小波系数的范围映射到0~255区间内
					short temp;
					temp = pTran1[iHstart + y][iWstart + x];
					temp -= MinPixVal;
					temp *= 255;
					if (Diff == 0)
						temp = 0;
					else
						temp /= Diff;
					spOriginData[iHstart + y][iWstart + x] = temp;

				}

			// 计算每个子频带的方差
			Vr[(i - 1) * 4] = ComputeVariance(iHstart, iHstart + iHeight_H,
					iWstart, iWstart + iWidth_H, pData);// LL子频带的方差
			Vr[(i - 1) * 4 + 1] = ComputeVariance(iHstart, iHstart + iHeight_H,
					iWstart + iWidth_H, iWstart + iWidth, pData);// HL子频带的方差
			Vr[(i - 1) * 4 + 2] = ComputeVariance(iHstart + iHeight_H, iHstart
					+ iHeight, iWstart, iWstart + iWidth_H, pData);// LH子频带的方差
			Vr[(i - 1) * 4 + 3] = ComputeVariance(iHstart + iHeight_H, iHstart
					+ iHeight, iWstart + iWidth_H, iWstart + iWidth, pData);// HH子频带的方差

			iHeight = iHeight >> 1;
			iWidth = iWidth >> 1;
			iHeight_H = iHeight / 2;
			iWidth_H = iWidth / 2;
		}

		return true;
	}

	// 计算矩阵p的方差,w为矩阵的宽，h为高
	public double ComputeVariance(int Hstart, int Hend, int Wstart, int Wend,
			short[][] LL) {
		double v = 0;
		double temp1 = 0;
		double temp2 = 0;
		double mean = 0;

		// 计算一阶矩
		for (int i = Hstart; i < Hend; i++)
			for (int j = Wstart; j < Wend; j++) {
				temp1 += (double) LL[i][j];
			}
		mean = temp1 / ((Hend - Hstart + 1) * (Wend - Wstart + 1));

		// 计算方差
		for (int i = Hstart; i < Hend; i++)
			for (int j = Wstart; j < Wend; j++) {
				temp2 += Math.pow(((double) LL[i][j] - mean), 2.0);

			}
		v = temp2 / ((Hend - Hstart + 1) * (Wend - Wstart + 1));

		return v;
	}

	public boolean Haar_Once(short[][] spOriginData, short[][] spTransData0,
			short[][] spTransData1, int Hstart, int Wstart, int nHeight,
			int nHeight_H, int nWidth, int nWidth_H, int layer) {
		
		int Trans_W, // 图像扫描线控制：横坐标
		Trans_H, // 图像扫描线控制：纵坐标
		Trans_M, // 图像矩阵的横坐标
		Trans_N; // 图像矩阵的纵坐标
		short Trans_Coeff0; // 小波变换系数
		short Trans_Coeff1; // 变换滤波系数
		// 本模块完成变换系数的赋值采样
		// 行变换,第一次（layer=1时）时nHeight即为原始图像的高度值

		if (layer == 1)// 如果是第一层，则数据源为spOriginData
		{

			for (Trans_H = 0; Trans_H < nHeight; Trans_H++) {

				for (Trans_N = 0; Trans_N < nWidth_H; Trans_N++) {
					Trans_W = Trans_N << 1;
					short temp1 = spOriginData[Trans_H + Hstart][Trans_W
							+ Wstart];
					short temp2 = spOriginData[Trans_H + Hstart][Trans_W
							+ Wstart + 1];
					// 计算平均值
					Trans_Coeff1 = (short) ((temp1 + temp2) >> 1);
					Trans_Coeff0 = (short) ((temp1 - temp2) >> 1);
					// 系数预测
					spOriginData[Trans_H + Hstart][Trans_W + Wstart] = Trans_Coeff1;
					spOriginData[Trans_H + Hstart][Trans_W + Wstart + 1] = Trans_Coeff0;
				}
				// 完成边界处理
				Trans_Coeff0 = (short) ((spOriginData[Trans_H + Hstart][Wstart + nWidth
						- 2] - spOriginData[Trans_H + Hstart][Wstart + nWidth
						- 1]) >> 1);
				spOriginData[Trans_H + Hstart][Wstart + nWidth - 1] = Trans_Coeff0;

				// 经过上述行处理后，偶数列为平均值，奇数列为细节值，因此将所有偶数列存储到前半部分，奇数列存储在后半部分
				for (Trans_N = 0; Trans_N < nWidth_H; Trans_N++) {
					Trans_W = Trans_N << 1;

					spTransData0[Trans_H + Hstart][Trans_N + Wstart] = (spOriginData[Trans_H
							+ Hstart][Trans_W + Wstart]);
					spTransData0[Trans_H + Hstart][nWidth_H + Trans_N + Wstart] = (spOriginData[Trans_H
							+ Hstart][Trans_W + 1 + Wstart]);

				}
				spTransData0[Trans_H + Hstart][Wstart + nWidth - 1] = spOriginData[Trans_H
						+ Hstart][Wstart + nWidth - 1];

			}
		}
		if (layer > 1)// 如果不是第一层，则数据源为spTransData1
			for (Trans_H = 0; Trans_H < nHeight; Trans_H++) {

				for (Trans_N = 0; Trans_N < nWidth_H; Trans_N++) {
					Trans_W = Trans_N << 1;
					short temp1 = spTransData1[Trans_H + Hstart][Trans_W
							+ Wstart];
					short temp2 = spTransData1[Trans_H + Hstart][Trans_W
							+ Wstart + 1];
					// 计算平均值
					Trans_Coeff1 = (short) ((temp1 + temp2) >> 1);
					Trans_Coeff0 = (short) ((temp1 - temp2) >> 1);
					// 系数预测
					spTransData1[Trans_H + Hstart][Trans_W + Wstart] = Trans_Coeff1;
					spTransData1[Trans_H + Hstart][Trans_W + Wstart + 1] = Trans_Coeff0;
				}
				// 完成边界处理
				Trans_Coeff0 = (short) ((spTransData1[Trans_H + Hstart][Wstart + nWidth
						- 2] - spTransData1[Trans_H + Hstart][Wstart + nWidth
						- 1]) >> 1);
				spTransData1[Trans_H + Hstart][Wstart + nWidth - 1] = Trans_Coeff0;

				// 经过上述行处理后，偶数列为平均值，奇数列为细节值，因此将所有偶数列存储到前半部分，奇数列存储在后半部分
				for (Trans_N = 0; Trans_N < nWidth_H; Trans_N++) {
					Trans_W = Trans_N << 1;

					spTransData0[Trans_H + Hstart][Trans_N + Wstart] = (spTransData1[Trans_H
							+ Hstart][Trans_W + Wstart]);
					spTransData0[Trans_H + Hstart][nWidth_H + Trans_N + Wstart] = (spTransData1[Trans_H
							+ Hstart][Trans_W + Wstart + 1]);

				}
				spTransData0[Trans_H + Hstart][Wstart + nWidth - 1] = spTransData1[Trans_H
						+ Hstart][Wstart + nWidth - 1];

			}

		// 水平方向的变换结束

		// 行提升后的数据在spTransData0中，spTransData0中的数据自然奇偶有序
		for (Trans_W = 0; Trans_W < nWidth; Trans_W++) {

			// 第一次提升奇数坐标系数
			for (Trans_M = 0; Trans_M < nHeight_H; Trans_M++) {
				// 竖直方向的变换
				Trans_H = Trans_M << 1;
				short temp1 = spTransData0[Trans_H + Hstart][Trans_W + Wstart];
				short temp2 = spTransData0[Trans_H + Hstart + 1][Trans_W
						+ Wstart];
				// 计算平均值
				Trans_Coeff1 = (short) ((temp1 + temp2) >> 1);
				Trans_Coeff0 = (short) ((temp1 - temp2) >> 1);
				// 系数预测
				spTransData0[Trans_H + Hstart][Trans_W + Wstart] = Trans_Coeff1;
				spTransData0[Trans_H + Hstart + 1][Trans_W + Wstart] = Trans_Coeff0;
			}
			// 完成边界处理
			Trans_Coeff0 = (short) ((spTransData0[Hstart + nHeight - 2][Wstart
					+ Trans_W] - spTransData1[Hstart + nHeight - 1][Wstart
					+ Trans_W]) >> 1);
			spTransData0[Hstart + nHeight - 1][Wstart + Trans_W] = Trans_Coeff0;

			// 经过上述列处理后，偶数行为平均值，奇数行为细节值，因此将所有偶数行存储到前半部分，奇数行存储在后半部分
			for (Trans_M = 0; Trans_M < nHeight_H; Trans_M++) {
				Trans_H = Trans_M << 1;

				spTransData1[Trans_M + Hstart][Trans_W + Wstart] = (spTransData0[Trans_H
						+ Hstart][Trans_W + Wstart]);
				spTransData1[nHeight_H + Trans_M + Hstart][Trans_W + Wstart] = (spTransData0[Trans_H
						+ Hstart + 1][Trans_W + Wstart]);

			}
			spTransData1[Hstart + nHeight - 1][Wstart + Trans_W] = spTransData0[Hstart
					+ nHeight - 1][Wstart + Trans_W];

		}

		return true;

	}

	

	//*******************************************************************************************//
	// Judge whether the frame is BlackWhiteFrame
	// Implemented by Tang Sheng 2007.7.25
	//*******************************************************************************************//
	//!此处函数未完成，原函数声明为 double VisFeaExtract::MeanStd(int []vector, int size, double  & std)
	public double MeanStd(int []vector, int size, double  std)
	{
		double Sum=0.0;
		for(int i=0;i<size;i++)
		{
			Sum=Sum+vector[i];
		}
		double mean=Sum/((double)(size));

		Sum=0.0;
		for(int i=0;i<size;i++)
		{
			Sum=Sum+(vector[i]-mean)*(vector[i]-mean);
		}
		Sum=Sum/(size-1);

		std = Math.sqrt(Sum);
		return mean;
	}
	//!此处函数未完成，原函数声明为 bool CVisFeaExtract::IsBlackWhiteFrame(double &meanU, double &stdU, double &meanV, double &stdV)
	public boolean IsBlackWhiteFrame(double meanU, double stdU, double meanV, double stdV)
	{
		if (m_pImage==null) 
			return false;

		int size=m_height*m_width;

		int []U=new int [size];
		int []V=new int [size];

		int count=0;
		CvMat mat = m_pImage.asCvMat();

		for(int i=0;i<m_height;i++)
		{
			for(int j=0;j<m_width;j++)
			{
				int R = (int) mat.get(i, j, 0);
				int G = (int) mat.get(i, j, 1);
				int B = (int) mat.get(i, j, 2);
				
				int Y= (int)(0.299f * R + 0.587f * G + 0.114f * B);
				U[count] = (int)((B-Y) * 0.565f + 128);
				V[count] = (int)((R-Y) * 0.713f + 128);
				count++;	
			}
		}

		meanU=MeanStd(U, size, stdU);
		meanV=MeanStd(V, size, stdV);

		if ( meanU>120 && meanU<136 &&
			 meanV>120 && meanV<136 &&
			 stdU + stdV < 5.0
			)
			return true;
		else
			return false;
	}

	
	
	
	private double max(double a, double b) {
		return a > b ? a : b;
	}

	private double min(double a, double b) {
		return a < b ? a : b;
	}

	public static void main(String[] args) throws IOException {
		// IplImage image = IplImage.createFrom(ImageIO
		// .read(new File("kebiao.png")));
		// CvMat mat = image.asCvMat();
		// int size = image.width() * image.height();
		// int nbChannels = mat.channels();
		// // Random random = new Random();
		// for (int i = 0; i < mat.rows(); i++) {
		// for (int j = 0; j < mat.cols(); j++)
		// for (int k = 0; k < nbChannels; k++)
		// mat.put(i, j, k, 255);
		// }
		// // System.out.println(mat.rows() + " "+ mat.cols() + " " +
		// image.width()
		// // + " " + image.height());
		// ImageIO.write(mat.asIplImage().getBufferedImage(), "png", new File(
		// "kebiao2.png"));

		double a[][][] = { { { 0 } } };
		System.out.println(a[0][0][0]);
		// test(a);
		System.out.println(a[0][0][0]);
	}

}
